home *** CD-ROM | disk | FTP | other *** search
/ CD ROM Paradise Collection 4 / CD ROM Paradise Collection 4 1995 Nov.iso / program / nrpas13.zip / RKQC.PAS < prev    next >
Pascal/Delphi Source File  |  1991-04-29  |  2KB  |  62 lines
  1. PROCEDURE rkqc(VAR y,dydx: glarray; n: integer; VAR x: real;
  2.       htry,eps: real; yscal: glarray; VAR hdid,hnext: real);
  3. (* Programs using routine RKQC must provide a
  4. PROCEDURE derivs(x:real; y:glnarray; VAR dydx:glnarray);
  5. which returns the derivatives dydx at location x, given both x and the
  6. function values y. They must also define the type
  7. TYPE
  8.    glarray = ARRAY [1..n] OF real;
  9. in the main routine.   *)
  10. LABEL 1;
  11. CONST
  12.    pgrow=-0.20;
  13.    pshrnk=-0.25;
  14.    fcor=0.06666666;   (* 1.0/15.0 *)
  15.    one=1.0;
  16.    safety=0.9;
  17.    errcon=6.0e-4;
  18. VAR
  19.    i: integer;
  20.    xsav,hh,h,temp,errmax: real;
  21.    dysav,ysav,ytemp: glarray;
  22. BEGIN
  23.    xsav := x;
  24.    FOR i := 1 TO n DO BEGIN
  25.       ysav[i] := y[i];
  26.       dysav[i] := dydx[i]
  27.    END;
  28.    h := htry;
  29. 1:   hh := 0.5*h;
  30.    rk4(ysav,dysav,n,xsav,hh,ytemp);
  31.    x := xsav+hh;
  32.    derivs(x,ytemp,dydx);
  33.    rk4(ytemp,dydx,n,x,hh,y);
  34.    x := xsav+h;
  35.    IF (x = xsav) THEN BEGIN
  36.       writeln('pause in routine RKQC');
  37.       writeln('stepsize too small'); readln
  38.    END;
  39.    rk4(ysav,dysav,n,xsav,h,ytemp);
  40.    errmax := 0.0;
  41.    FOR i := 1 TO n DO BEGIN
  42.       ytemp[i] := y[i]-ytemp[i];
  43.       temp := abs(ytemp[i]/yscal[i]);
  44.       IF (errmax < temp) THEN errmax := temp
  45.    END;
  46.    errmax := errmax/eps;
  47.    IF (errmax > one) THEN BEGIN
  48.       h := safety*h*exp(pshrnk*ln(errmax));
  49.       GOTO 1 END
  50.    ELSE BEGIN
  51.       hdid := h;
  52.       IF (errmax > errcon) THEN BEGIN
  53.          hnext := safety*h*exp(pgrow*ln(errmax))
  54.       END ELSE BEGIN
  55.          hnext := 4.0*h
  56.       END
  57.    END;
  58.    FOR i := 1 TO n DO BEGIN
  59.       y[i] := y[i]+ytemp[i]*fcor
  60.    END
  61. END;
  62.